;--------------------------------------------------------------------------------------------------
; Arquivo: gmm.pro (funcao)
; Objetivo:
;   
; Entradas: 
;  
; Saidas:
;  
; Observacoes:
;   
; Autor: Daniela Souza dos Anjos
; Data: 08/09/2011
; Historico de Mudanças:
; 
;--------------------------------------------------------------------------------------------------PRO gmm


pro teste
vMi=[0,0,0]
vPadrao=[1,1,1]
mCovInicial=NovoSigma(2,3)
print,ProbDeXSerDaClasse2(vMi,vPadrao,mCovInicial[*,*,0])

end


pro executar, nIter, nclasses
imagem=carregarImagem()
outfile = DIALOG_PICKFILE(FILTER='*.tif',/WRITE)
pixels=imagem.pixelsImagem
npixels=n_elements(pixels[*,0])
nbandas=n_elements(pixels[0,*])
limiteMax=max(pixels)
NovaMedia=fix(randomu(seed,nclasses,nbandas)*limiteMax)
mNovoSigma=NovoSigma(nclasses,nbandas)*(limiteMax/12.)^2.
;mNovoSigma[*]=255./24.
print, "Iniciando as estimativas"
Resultado=Calcgmm(pixels,mNovoSigma,NovaMedia,nIter)
print, "Classificacao MAP"
imgClassificada=map(pixels,Resultado.mMi,Resultado.mSigma, Resultado.vProb)
print, imgClassificada[1:4,*]


saida=reform(imgClassificada,imagem.ncolunas,npixels/imagem.ncolunas)
WRITE_TIFF,outfile,saida,/SHORT
print, "Finalizado."
end



function MAP, mPadrao, mMi, mSigma, vProb
  npixels=n_elements(mPadrao[*,0])
  nbandas=n_elements(mPadrao[0,*])
  nclasses=n_elements(vProb)
  imgClassificada=intarr(npixels)
  ;ProbAPosteriori=fltarr(npixels,nclasses)
  for i=0L, npixels-1 do begin
    classeMax=0
    PjxMax=0
    for j=0, nclasses-1 do begin
      Pjx=vProb[j]*ProbDeXSerDaClasse(mMi[j,*],mPadrao[i,*], mSigma[*,*,j])
      if Pjx gt PjxMax then begin
        PjxMax=Pjx
        classeMax=j
      end
    end
    imgClassificada[i]=classeMax
  end
  return, imgClassificada
end

function ProbDeXSerDaClasse, vMi, vPadrao, mCov 
;x dado j (14.12 T&K)
;Calcula a probabilidade do padrao dado pertencer a distribuicao com media vMi e covariancia mCov
;vMi - vetor de media da classe para cada banda
;vPadrao - vetor com os valores de pixel em cada banda
;mCov - matriz de covariancia
nbandas=n_elements(vMi)
expoente= (-0.5*(vPadrao-vMi) # invert(mCov) # transpose(vPadrao-vMi))[0]
denominador=sqrt((2.*!pi)^(nbandas) * norm(mCov))
P=exp( expoente)/denominador
if ~finite(P) then begin
  p=0
end
return, P
end

;function ProbDeXSerDaClasse, vMi, vPadrao, fSigma 
;;x dado j (2.96 T&K)
;;Calcula a probabilidade do padrao dado pertencer a distribuicao com media vMi e covariancia Sigma*I
;;vMi - vetor de media da classe para cada banda
;;vPadrao - vetor com os valores de pixel em cada banda
;;fSigma - float - variancia da classe
;nbandas=n_elements(vMi)
;expoente=exp(-((norm(vPadrao-vMi))^2.)/(2.*fSigma^2.))
;P=(((2.*!pi*fsigma^2.)^(-nbandas/2.)))* expoente
;if ~finite(P) then begin
;  p=0
;end
;return, P
;end

function ProbDaClasseSerDeX, mMi, mSigma, idJ, vPadrao,vProb
;j dado x (Bayes - 2.101 T&K)
;mMI - matriz com os vetores de medias das classes para cada banda (n classes X nbandas)
;vSigma - vetor com as variancias de cada classe (1 X nclasses)
;vPadrao - vetor com os valores de pixel em cada banda (1 X npixels)
;vProb - vetor com as probabilidades de cada classe (1 X nclasses)
pjx=ProbDeXSerDaClasse(mMi[idJ,*], vPadrao, mSigma[*,*,idJ])*vProb[idJ]
px=ProbDEX(mMi, vPadrao, mSigma,vProb)
P=pjx/px
if ~finite(P) then begin
  P=0.
end
return, p
end

function ProbDeX,mMi, vPadrao, mSigma,vProb
;2.102 T&K
j=n_elements(vProb)
p=0.
for i=0,j-1 do begin
p+=ProbDeXSerDaClasse(mMi[i,*], vPadrao, mSigma[*,*,i])*vProb[i]
endfor
return, p
end

function Maximization, mMI, mSigma, mPadrao, vProb
;2.98 T&K/;2.99 T&K/;2.100 T&K
npixels=n_elements(mPadrao[*,0])
nbandas=n_elements(mPadrao[0,*])
nclasses=n_elements(vProb)
mNovoSigma=fltarr(nbandas,nbandas,nclasses)
mNovomi=fltarr(nclasses,nbandas)
vNovoProb=fltarr(nclasses)

for idJ=0,nclasses-1 do begin
  vp1=transpose(dblarr(nbandas))
  p2=double(0)
  mS1=dblarr(nbandas,nbandas)
  for k=0L,npixels-1 do begin
    ; estima a media
    Pjx=ProbDaClasseSerDeX( mMi, mSigma, idJ, mPadrao[k,*],vProb)
    ;if (k mod 1000) eq 0 then begin
    if (( (k+idj*npixels)*100./(npixels*nclasses)) mod 1) eq 0 then begin
      print, (k+idj*npixels)*100/(npixels*nclasses)
    end
    dist=mPadrao[k,*]-mMi[idJ,*]
    mS1+=Pjx*(transpose(dist) # dist)
    vp1+=Pjx* mPadrao[k,*]
    p2+=Pjx
  endfor

  if p2 ne 0 then begin
    ;estima novos parametros apenas se a probabilidade somada de todos os padroes for diferente de zero.
    mNovoMi[idj,*]=vp1/p2
    vNovoProb[idj]=p2/npixels
    mNovoSigma[*,*,idj]=mS1/p2
  end else begin
    ;se a probabilidade eh igual a zero, refazer a estimacao da media inicial.
    limite=max(mPadrao)
    mNovoMi[idj,*]=fix(randomu(seed,nbandas)*limite)
    vNovoProb[idj]=vProb[idj]
    mNovoSigma[*,*,idj]=limite/12.
  end
  
    
endfor


resultado={mMi:mNovoMi,mSigma:mNovoSigma, vProb:vNovoProb}
return, resultado
end  



function NovoSigma, NumeroClasses, NumeroAtributos
  ;sigma0=[]
  sigma0=identity(NumeroAtributos)
  for i=1, NumeroClasses-1 do begin
    sigma0=[[sigma0],[identity(NumeroAtributos)]]
  endfor
  
  return, reform(sigma0,NumeroAtributos,NumeroAtributos,NumeroClasses)
end
 
  
function CalcGmm, mPadrao, SigmaInicial, MediaInicial, niter
;ProbInicial
;Dados de entrada
;Numero de classes(inteiro)
;SigmaInicial - Vetor de matrizes de covariancia, uma matriz pra cada classes(Matriz 3D - classe i, atributo j, atributo k)
;MediaInicial - Media de cada classe em cada banda (Matriz 2D - classe i, banda n)
;ProbInicial - Probabilidade a priori de cada classe (Vetor com um num real pra cada classe)
nclasses=n_elements(MediaInicial[*,0])
ProbInicial=MAKE_ARRAY(nclasses,/FLOAT, VALUE =1./nclasses)
iteracao={mMi:MediaInicial,mSigma:SigmaInicial,vProb:ProbInicial}
for i=0, niter-1 do begin
print, "Iteracao "+strtrim(i+1,1)
Iteracao=Maximization( iteracao.mMi, Iteracao.mSigma, mPadrao, Iteracao.vProb)
endfor
return, Iteracao
;iteracao1.mMI
end



pro gmm
;iteracoes, classes
executar, 1, 5
;teste
END